Research Questions:

  1. How does stream metabolism vary based on catchment characteristics, antecedent flow conditions, and nitrogen availability?

    • Figures 2, 3
  2. Is variance in nitrogen demand and instream uptake (either as ammonium or nitrate) more related to biological (GPP, ER, biomass or chlorophyll-a), physical processes (flow and water temperature), or ambient nutrient and organic matter availability?

    • Figures 4, 5
  3. How does in-stream productivity and nitrogen supply and demand dynamics change within stream catchments and across dry (2021 and 2022) and wet water years (2023)?

    • Figures 6

Figure 1: Map

Not included here.

Pause to make a column for biomass to chla ratios And also make biomass in in micro grams to more comparable to chl-a

Question 1

How does stream metabolism vary based on catchment characteristics, antecedent flow conditions, and nitrogen availability?

Figure 2: Metabolism Time series

  • The magnitude of GPP and ER across all sites

Where GPP is in green, ER is in purple, NEP is in black dataframe: metab_dat

Time series metabolism plots for each site

Descriptive stats to describe stream based differences in

Both BW to GB - large to small streams, as well as variation with stream.

  • Water temp

  • Flow

  • OM

  • Biomass

  • Chla

  • GPP

  • |ER|

BW catchment:

##   mean_AFDM_mgg min_AFDM_mgg max_AFDM_mgg mean_Chla_ugL_Q min_Chla_ugL_Q
## 1        16.642     8.567217     49.13182          45.981      0.3621839
##   max_Chla_ugL_Q mean_AFDM_ugcm2 min_AFDM_ugcm2 max_AFDM_ugcm2 mean_GPP min_GPP
## 1       417.6117          51.968              0            369    1.883   0.001
##    max_GPP mean_ER    min_ER    max_ER mean_wtr_m min_wtr_m max_wtr_m mean_Q_m
## 1 11.54822 -12.379 -28.72293 -1.725866      8.614 0.2916944   20.4719    0.715
##       min_Q_m  max_Q_m
## 1 0.009937639 8.531119
## # A tibble: 2 × 16
##   substrate mean_PO4_ugL_dl min_PO4_ugL_dl max_PO4_ugL_dl mean_NO3_mgL_dl
##   <chr>               <dbl>          <dbl>          <dbl>           <dbl>
## 1 pw                  10.5            4.26           62.3           0.063
## 2 sw                   9.65           1.68           26.2           0.015
## # ℹ 11 more variables: min_NO3_mgL_dl <dbl>, max_NO3_mgL_dl <dbl>,
## #   mean_NH4_mgL_dl <dbl>, min_NH4_mgL_dl <dbl>, max_NH4_mgL_dl <dbl>,
## #   mean_pH <dbl>, min_pH <dbl>, max_pH <dbl>, mean_DOC_mgL_dl <dbl>,
## #   min_DOC_mgL_dl <dbl>, max_DOC_mgL_dl <dbl>

GB Catchment

##   mean_AFDM_mgg min_AFDM_mgg max_AFDM_mgg mean_Chla_ugL_Q min_Chla_ugL_Q
## 1        22.175     8.710089     58.84203          48.077       2.046401
##   max_Chla_ugL_Q mean_AFDM_ugcm2 min_AFDM_ugcm2 max_AFDM_ugcm2 mean_GPP
## 1       260.7374          86.347              0            292    0.107
##        min_GPP  max_GPP mean_ER    min_ER     max_ER mean_wtr_m  min_wtr_m
## 1 4.520732e-05 3.150141   -4.97 -19.23026 -0.4027769      7.826 0.01800694
##   max_wtr_m mean_Q_m      min_Q_m   max_Q_m
## 1  15.66177    0.054 0.0009168555 0.9424983
## # A tibble: 2 × 16
##   substrate mean_PO4_ugL_dl min_PO4_ugL_dl max_PO4_ugL_dl mean_NO3_mgL_dl
##   <chr>               <dbl>          <dbl>          <dbl>           <dbl>
## 1 pw                   26.7           8.5           109.            0.084
## 2 sw                   14.1           8.47           25.3           0.017
## # ℹ 11 more variables: min_NO3_mgL_dl <dbl>, max_NO3_mgL_dl <dbl>,
## #   mean_NH4_mgL_dl <dbl>, min_NH4_mgL_dl <dbl>, max_NH4_mgL_dl <dbl>,
## #   mean_pH <dbl>, min_pH <dbl>, max_pH <dbl>, mean_DOC_mgL_dl <dbl>,
## #   min_DOC_mgL_dl <dbl>, max_DOC_mgL_dl <dbl>

Break out by stream reach location:

BWL:
BWU:
GBL:
GBU:
Kernal density plots of GPP and ER

Colored by years

Modified figure 2

(Q1) - How does stream metabolism vary based on catchment characteristics, antecedent flow conditions, and nitrogen availability?

General drivers of productivity

Both standing stocks

  • A. Epilithic biomass

  • B. Epilithic chl-a

Ecosystem fluxes

  • C. Gross primary productivity (GPP)

  • E. Absolute value of ecosystem respiration (ER)

Each productivity response is individually modeled within stream using glm

dataframe = covariat_datq

Fixed effect co-variance ?

Maybe some with flow and water temp but all correlations are under 0.6

At Blackwood:

cor_data <- covariat_datq_BW[, c(6, 11:14, 18, 28, 30)]

chart.Correlation(cor_data, histogram=TRUE, pch=19)

At Glenbrook:

cor_data <- covariat_datq_GB[, c(6, 11:14, 18, 28, 30)]

chart.Correlation(cor_data, histogram=TRUE, pch=19)

A. Biomass models

BW:

biomass_mod_bw <- lmer(log(AFDM_ugcm2+1)~ 
                         scale(NO3_mgL_dl* 1000) + 
                         scale(NH4_mgL_dl* 1000) + 
                         scale(PO4_ugL_dl) + 
                        scale(AFDM_mgg)+
                        # scale(DOC_mgL_dl)+
                         scale(wtr_m) + 
                         scale(Q_m) + (1|site), data=covariat_datq_BW)

summary(biomass_mod_bw)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(AFDM_ugcm2 + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_BW
## 
## REML criterion at convergence: 231.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.14121 -0.25057  0.02282  0.43491  2.21972 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.07452  0.273   
##  Residual             1.81282  1.346   
## Number of obs: 66, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)               3.46061    0.29677  0.97047  11.661  0.05825 . 
## scale(NO3_mgL_dl * 1000)  0.23732    0.15059 58.86194   1.576  0.12040   
## scale(NH4_mgL_dl * 1000) -0.19969    0.18151 58.07195  -1.100  0.27581   
## scale(PO4_ugL_dl)         0.08077    0.12331 58.59587   0.655  0.51500   
## scale(AFDM_mgg)           0.18559    0.29544 58.06074   0.628  0.53235   
## scale(wtr_m)              0.68938    0.25854 58.00248   2.666  0.00992 **
## scale(Q_m)               -0.07222    0.15056 58.94510  -0.480  0.63326   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10  0.022                                   
## s(NH4_L_*10 -0.119 -0.356                            
## scl(PO4_L_)  0.012 -0.131 -0.045                     
## scl(AFDM_m) -0.068 -0.024 -0.071  0.177              
## scal(wtr_m)  0.290  0.229  0.074 -0.027 -0.169       
## scale(Q_m)  -0.159  0.196  0.243 -0.253 -0.201  0.410
vif(biomass_mod_bw)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 1.321235                 1.298738                 1.112896 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.076703                 1.269469                 1.435108
#hist(residuals(biomass_mod_bw), main = paste(NULL))
r.squaredGLMM(biomass_mod_bw)
##            R2m       R2c
## [1,] 0.1554381 0.1887849

GB:

biomass_mod_gb <- lmer(log(AFDM_ugcm2+1)~ 
                         scale(NO3_mgL_dl* 1000) + 
                         scale(NH4_mgL_dl* 1000) + 
                         scale(PO4_ugL_dl) + 
                         scale(AFDM_mgg)+
                         #scale(DOC_mgL_dl)+
                         scale(wtr_m) + 
                         scale(Q_m)+(1|site), data=covariat_datq_GB)

summary(biomass_mod_gb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(AFDM_ugcm2 + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_GB
## 
## REML criterion at convergence: 207.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95119 -0.45345  0.01751  0.57955  1.76719 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.000    0.000   
##  Residual             1.592    1.262   
## Number of obs: 61, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)               3.70075    0.20246 54.00000  18.279   <2e-16 ***
## scale(NO3_mgL_dl * 1000) -0.61361    0.29306 54.00000  -2.094   0.0410 *  
## scale(NH4_mgL_dl * 1000)  0.28855    0.13960 54.00000   2.067   0.0435 *  
## scale(PO4_ugL_dl)        -0.03695    0.12569 54.00000  -0.294   0.7699    
## scale(AFDM_mgg)          -0.38478    0.18041 54.00000  -2.133   0.0375 *  
## scale(wtr_m)              0.04233    0.22999 54.00000   0.184   0.8547    
## scale(Q_m)                0.02593    0.13086 54.00000   0.198   0.8437    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.123                                   
## s(NH4_L_*10 -0.157 -0.302                            
## scl(PO4_L_) -0.082 -0.035 -0.326                     
## scl(AFDM_m) -0.275  0.146 -0.310  0.055              
## scal(wtr_m) -0.314  0.354  0.020 -0.100 -0.124       
## scale(Q_m)  -0.441  0.316  0.129 -0.141  0.165  0.387
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
vif(biomass_mod_gb)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 1.399768                 1.468211                 1.150803 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.230784                 1.346312                 1.393157
#hist(residuals(biomass_mod_gb), main = paste(NULL))
r.squaredGLMM(biomass_mod_gb)
##           R2m      R2c
## [1,] 0.139549 0.139549

B. Chlorophyll-a models

BW:

chla_mod_bw <- lmer(log(Chla_ugL_Q+1)~ 
                         scale(NO3_mgL_dl* 1000) + 
                         scale(NH4_mgL_dl* 1000) + 
                         scale(PO4_ugL_dl)+ 
                     scale(AFDM_mgg)+
                         #scale(DOC_mgL_dl)+
                         scale(wtr_m) + 
                         scale(Q_m)+(1|site), data=covariat_datq_BW)

summary(chla_mod_bw)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Chla_ugL_Q + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_BW
## 
## REML criterion at convergence: 256.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7932 -0.5506 -0.1391  0.4846  2.4774 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.4873   0.6981  
##  Residual             2.2105   1.4868  
## Number of obs: 69, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)               2.69172    0.54445  1.05614   4.944   0.1173  
## scale(NO3_mgL_dl * 1000)  0.26748    0.18377 61.29299   1.456   0.1506  
## scale(NH4_mgL_dl * 1000) -0.30309    0.16752 61.69467  -1.809   0.0753 .
## scale(PO4_ugL_dl)         0.14788    0.13070 61.25780   1.131   0.2623  
## scale(AFDM_mgg)           0.20959    0.29054 61.36433   0.721   0.4734  
## scale(wtr_m)              0.00381    0.31804 61.06972   0.012   0.9905  
## scale(Q_m)               -0.05790    0.18815 61.88681  -0.308   0.7593  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.132                                   
## s(NH4_L_*10  0.035 -0.360                            
## scl(PO4_L_)  0.066 -0.163 -0.052                     
## scl(AFDM_m)  0.189 -0.262  0.076  0.118              
## scal(wtr_m) -0.158  0.525 -0.220 -0.025 -0.228       
## scale(Q_m)  -0.199  0.434  0.021 -0.242 -0.370  0.618
vif(chla_mod_bw)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 1.655305                 1.245731                 1.115461 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.182510                 1.998142                 2.057583
#hist(residuals(chla_mod_bw), main = paste(NULL))
r.squaredGLMM(chla_mod_bw)
##            R2m       R2c
## [1,] 0.0826535 0.2483622

GB:

chla_mod_gb <- lmer(log(Chla_ugL_Q+1)~ 
                      scale(NO3_mgL_dl* 1000) + 
                      scale(NH4_mgL_dl* 1000) + 
                      scale(PO4_ugL_dl)+
                     scale(AFDM_mgg)+
                      #scale(DOC_mgL_dl)+
                      scale(wtr_m) + 
                      scale(Q_m) + (1|site), data=covariat_datq_GB)

summary(chla_mod_gb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Chla_ugL_Q + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_GB
## 
## REML criterion at convergence: 220.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9286 -0.7348  0.1622  0.4396  2.1911 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.02052  0.1432  
##  Residual             0.72393  0.8508  
## Number of obs: 83, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)               3.44327    0.16017  1.08780  21.497  0.02305 * 
## scale(NO3_mgL_dl * 1000) -0.35575    0.18547 75.95345  -1.918  0.05885 . 
## scale(NH4_mgL_dl * 1000)  0.50245    0.15957 72.39039   3.149  0.00238 **
## scale(PO4_ugL_dl)        -0.06480    0.08634 70.76286  -0.751  0.45538   
## scale(AFDM_mgg)          -0.15622    0.12522 39.22083  -1.248  0.21960   
## scale(wtr_m)              0.34064    0.14431 64.56896   2.361  0.02128 * 
## scale(Q_m)                0.17167    0.08714 75.96634   1.970  0.05248 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.238                                   
## s(NH4_L_*10  0.050 -0.518                            
## scl(PO4_L_) -0.091  0.035 -0.424                     
## scl(AFDM_m)  0.060  0.324 -0.519  0.150              
## scal(wtr_m) -0.450  0.465 -0.038 -0.147  0.006       
## scale(Q_m)  -0.317  0.355  0.008 -0.115  0.086  0.439
vif(chla_mod_gb)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 2.059524                 2.264356                 1.306441 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.406260                 1.514430                 1.336345
#hist(residuals(chla_mod_gb), main = paste(NULL))
r.squaredGLMM(chla_mod_gb)
##            R2m       R2c
## [1,] 0.2204337 0.2419218

C. GPP models

BW:

gpp_mod_bw <- lmer(log(GPP_mean+1) ~ 
                      scale(NO3_mgL_dl* 1000) + 
                      scale(NH4_mgL_dl* 1000) + 
                      scale(PO4_ugL_dl)+
                      scale(AFDM_mgg)+
                      #scale(DOC_mgL_dl)+
                      scale(wtr_m) +
                      scale(Q_m)+ (1|site), 
                  data=covariat_datq_BW_na)

summary(gpp_mod_bw)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(GPP_mean + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_BW_na
## 
## REML criterion at convergence: 74.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3520 -1.0599  0.1177  0.6725  2.2217 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.2085   0.4567  
## Number of obs: 47, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)               0.60126    0.07173 40.00000   8.382 2.41e-10 ***
## scale(NO3_mgL_dl * 1000) -0.03634    0.14040 40.00000  -0.259 0.797102    
## scale(NH4_mgL_dl * 1000)  0.03251    0.10707 40.00000   0.304 0.762951    
## scale(PO4_ugL_dl)         0.05057    0.07937 40.00000   0.637 0.527619    
## scale(AFDM_mgg)           0.12434    0.08933 40.00000   1.392 0.171664    
## scale(wtr_m)              0.47138    0.16981 40.00000   2.776 0.008331 ** 
## scale(Q_m)                0.41802    0.11104 40.00000   3.765 0.000536 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.012                                   
## s(NH4_L_*10 -0.089 -0.791                            
## scl(PO4_L_)  0.128 -0.527  0.353                     
## scl(AFDM_m)  0.260  0.122 -0.269  0.086              
## scal(wtr_m) -0.195  0.605 -0.477 -0.318 -0.028       
## scale(Q_m)  -0.105  0.125 -0.035 -0.013 -0.096  0.613
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
vif(gpp_mod_bw)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 3.943474                 2.987120                 1.453622 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.151354                 3.006878                 1.951788
#hist(residuals(gpp_mod_bw),main = paste(NULL))
r.squaredGLMM(gpp_mod_bw)
##            R2m       R2c
## [1,] 0.3142777 0.3142777

GB:

Highly skewed even when log transformed.

gpp_mod_gb <- lmer(log(GPP_mean+1) ~ 
                      scale(NO3_mgL_dl* 1000) + 
                      scale(NH4_mgL_dl* 1000) + 
                      scale(PO4_ugL_dl)+ 
                             scale(AFDM_mgg)+
                      #scale(DOC_mgL_dl)+
                      scale(wtr_m) + 
                      scale(Q_m) + (1|site), data=covariat_datq_GB_na)

summary(gpp_mod_gb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(GPP_mean + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_GB_na
## 
## REML criterion at convergence: -13.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2910 -0.6469 -0.1149  0.2579  2.9371 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  site     (Intercept) 0.0002591 0.0161  
##  Residual             0.0128017 0.1131  
## Number of obs: 31, groups:  site, 2
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)               0.093556   0.033210  0.489448   2.817   0.3841  
## scale(NO3_mgL_dl * 1000)  0.007607   0.025168 23.999945   0.302   0.7651  
## scale(NH4_mgL_dl * 1000)  0.032072   0.020539  3.214033   1.562   0.2103  
## scale(PO4_ugL_dl)        -0.056735   0.029337  8.430856  -1.934   0.0873 .
## scale(AFDM_mgg)          -0.003998   0.019357  8.228514  -0.207   0.8414  
## scale(wtr_m)              0.009400   0.036439 23.006449   0.258   0.7987  
## scale(Q_m)                0.004899   0.018976 23.055864   0.258   0.7986  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.044                                   
## s(NH4_L_*10  0.141 -0.054                            
## scl(PO4_L_) -0.295 -0.549 -0.336                     
## scl(AFDM_m) -0.382  0.257 -0.289  0.045              
## scal(wtr_m) -0.518  0.147 -0.079 -0.146  0.376       
## scale(Q_m)  -0.431  0.021 -0.067  0.133  0.220  0.334
vif(gpp_mod_gb)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 1.721912                 1.307889                 1.935489 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.378244                 1.340005                 1.190482
#hist(residuals(gpp_mod_gb), main = paste(NULL))
r.squaredGLMM(gpp_mod_gb)
##            R2m       R2c
## [1,] 0.1675341 0.1840497

fixed_effects_nep1 <- broom.mixed::tidy(gpp_mod_gb, effects = "fixed", conf.int = TRUE)
# Filter to exclude intercept and arrange predictors
fixed_effects_nep1 <- fixed_effects_nep1[fixed_effects_nep1$term != "(Intercept)", ]
fixed_effects_nep1 <- fixed_effects_nep1[order(fixed_effects_nep1$estimate), ]  # Order by estimate
fixed_effects_nep1$Catchment <- "GB"

D. |ER| models

BW:

er_mod_bw <- lmer(log(ab_ER+1) ~ 
                      scale(NO3_mgL_dl* 1000) + 
                      scale(NH4_mgL_dl* 1000) + 
                      scale(PO4_ugL_dl)+
                      scale(AFDM_mgg)+
                      #scale(DOC_mgL_dl)+
                      scale(wtr_m)  +
                      scale(Q_m) + (1|site), 
                  data=covariat_datq_BW_na)

summary(er_mod_bw)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(ab_ER + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_BW_na
## 
## REML criterion at convergence: -60.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.24094 -0.59081  0.01924  0.68879  3.15168 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.121445 0.34849 
##  Residual             0.006303 0.07939 
## Number of obs: 47, groups:  site, 2
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               2.565975   0.246736  0.995205  10.400 0.061662 .  
## scale(NO3_mgL_dl * 1000)  0.014010   0.024947 39.016569   0.562 0.577598    
## scale(NH4_mgL_dl * 1000) -0.012010   0.018634 39.000887  -0.645 0.523010    
## scale(PO4_ugL_dl)        -0.004434   0.015607 39.083279  -0.284 0.777816    
## scale(AFDM_mgg)           0.002081   0.016103 39.027033   0.129 0.897859    
## scale(wtr_m)              0.133738   0.034606 39.103309   3.865 0.000408 ***
## scale(Q_m)               -0.133371   0.020169 39.032426  -6.613  7.3e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.001                                   
## s(NH4_L_*10 -0.004 -0.783                            
## scl(PO4_L_)  0.004 -0.359  0.289                     
## scl(AFDM_m)  0.014  0.060 -0.247 -0.050              
## scal(wtr_m) -0.010  0.613 -0.432  0.004 -0.161       
## scale(Q_m)  -0.006  0.177 -0.047  0.125 -0.165  0.652
vif(er_mod_bw)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 3.994356                 2.919201                 1.267362 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.187484                 3.437483                 2.045603
#hist(residuals(er_mod_bw), main = paste(NULL))
r.squaredGLMM(er_mod_bw)
##            R2m       R2c
## [1,] 0.2068482 0.9608681

# Extract fixed effects and confidence intervals
fixed_effects_er <- broom.mixed::tidy(er_mod_bw, effects = "fixed", conf.int = TRUE)
# Filter to exclude intercept and arrange predictors
fixed_effects_er <- fixed_effects_er[fixed_effects_er$term != "(Intercept)", ]
fixed_effects_er <- fixed_effects_er[order(fixed_effects_er$estimate), ]  # Order by estimate
fixed_effects_er$Catchment <- "BW"

GB:

er_mod_gb <- lmer(log(ab_ER+1) ~ 
                      scale(NO3_mgL_dl* 1000) + 
                      scale(NH4_mgL_dl* 1000) + 
                      scale(PO4_ugL_dl)+
                      scale(AFDM_mgg)+
                     # scale(DOC_mgL_dl)+
                      scale(wtr_m) +
                      scale(Q_m)+ (1|site), 
                  data=covariat_datq_GB_na)

summary(er_mod_gb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(ab_ER + 1) ~ scale(NO3_mgL_dl * 1000) + scale(NH4_mgL_dl *  
##     1000) + scale(PO4_ugL_dl) + scale(AFDM_mgg) + scale(wtr_m) +  
##     scale(Q_m) + (1 | site)
##    Data: covariat_datq_GB_na
## 
## REML criterion at convergence: 30.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.97321 -0.72982  0.02548  0.61764  2.00362 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.66140  0.8133  
##  Residual             0.06701  0.2589  
## Number of obs: 31, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)               1.71037    0.58126  1.01003   2.943  0.20656    
## scale(NO3_mgL_dl * 1000) -0.04198    0.05876 23.03354  -0.715  0.48210    
## scale(NH4_mgL_dl * 1000)  0.11777    0.05831 23.27930   2.020  0.05507 .  
## scale(PO4_ugL_dl)        -0.07954    0.07675 23.19187  -1.036  0.31074    
## scale(AFDM_mgg)           0.04697    0.05074 23.19413   0.926  0.36414    
## scale(wtr_m)              0.27939    0.08671 23.06356   3.222  0.00376 ** 
## scale(Q_m)                0.28472    0.04344 23.00103   6.554 1.09e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s(NO*1 s(NH*1 s(PO4_ s(AFDM scl(_)
## s(NO3_L_*10 -0.021                                   
## s(NH4_L_*10  0.061 -0.160                            
## scl(PO4_L_) -0.072 -0.374 -0.524                     
## scl(AFDM_m) -0.082  0.317 -0.492  0.271              
## scal(wtr_m) -0.087  0.193 -0.224  0.011  0.450       
## scale(Q_m)  -0.054  0.014 -0.033  0.099  0.175  0.312
vif(er_mod_gb)
## scale(NO3_mgL_dl * 1000) scale(NH4_mgL_dl * 1000)        scale(PO4_ugL_dl) 
##                 1.655592                 1.912623                 2.126215 
##          scale(AFDM_mgg)             scale(wtr_m)               scale(Q_m) 
##                 1.773162                 1.395697                 1.137209
#hist(residuals(er_mod_gb), main = paste(NULL))
r.squaredGLMM(er_mod_gb)
##            R2m       R2c
## [1,] 0.1538996 0.9221604

fixed_effects_er1 <- broom.mixed::tidy(er_mod_gb, effects = "fixed", conf.int = TRUE)
# Filter to exclude intercept and arrange predictors
fixed_effects_er1 <- fixed_effects_er1[fixed_effects_er1$term != "(Intercept)", ]
fixed_effects_er1 <- fixed_effects_er1[order(fixed_effects_er1$estimate), ]  # Order by estimate
fixed_effects_er1$Catchment <- "GB"

Figure 3: General drivers for biomass, chl-a and GPP

(Q1) - How does stream metabolism vary based on catchment characteristics, antecedent flow conditions, and nitrogen availability?

Mainly addresses how productivity is associated with nutrients and water quality:

In general GPP may be more a function of physical conditions but is partially associated with NH4.

Part2/ Question 2

Nitrogen demand and supply dynamics

Calculate nitrogen supply and demand:

# Constants
ra <- 0.5  # Autotrophic respiration coefficient (Hall & Tank 2003)
C_Nauto <- 16  # Autotrophic C:N ratio (Stelzer & Lamberti 2001)
C_Nhetero <- 20  # Heterotrophic C:N ratio (Hall & Tank 2003)
HGE <- 0.05  # Heterotrophic Growth Efficiency (Hall & Tank 2003)

# Calculate components of nitrogen demand
covariat_nitrogen2 <- covariat_datq %>%
  filter(substrate=="sw")%>% # select just surface water samples
  group_by(site) %>% 
  arrange(date) %>%  # Ensure chronological order
  mutate( ## 17 NAs
    v_mi = ifelse(is.na(v_m), 
                  rollapplyr(v_m, width = 3, FUN = function(x) mean(x, na.rm = TRUE), fill = NA, partial = TRUE), 
                  v_m),
    w_mi = ifelse(is.na(w_m), 
                  rollapplyr(w_m, width = 3, FUN = function(x) mean(x, na.rm = TRUE), fill = NA, partial = TRUE), 
                  w_m)
  ) %>%
  # Carry forward any remaining NAs
  mutate(
    v_mi = zoo::na.locf(v_mi, na.rm = FALSE),
    v_mi = zoo::na.locf(v_mi, fromLast = TRUE, na.rm = FALSE),
    w_mi = zoo::na.locf(w_mi, na.rm = FALSE),
    w_mi = zoo::na.locf(w_mi, fromLast = TRUE, na.rm = FALSE)
  ) %>%
  # If there are still missing values, replace them with the site mean
  mutate(
    v_mi = ifelse(is.na(v_mi), mean(v_m, na.rm = TRUE), v_mi),
    w_mi = ifelse(is.na(w_mi), mean(w_m, na.rm = TRUE), w_mi)
  ) %>%
  ungroup()%>%
  ## should already be at detection limits
  # filter(water_year < 2024) %>%
  # mutate(NO3_mgL_dl = if_else(NO3_mgL_dl < 0.0015, 0.0015, NO3_mgL_dl))%>%
  #   mutate(NH4_mgL_dl=if_else(NH4_mgL_dl < 0.001, 0.001, NH4_mgL_dl))%>%
  #   mutate(PO4_ugL_dl=if_else(PO4_ugL_dl < 0.201, 0.201, PO4_ugL_dl))%>%
  mutate(
    
    ## ** Pause here to double check Q units
    Q_Ls = Q_m * 1000, # flow from csm to Ls
    # Autotrophic respiration (raGPP)
    raGPP = GPP_mean * ra,
    # Autotrophic assimilation of N
    Auto_N_assim = GPP_mean / C_Nauto,
    # Heterotrophic respiration (Rh)
    Rh = ER_mean - raGPP,
    # Heterotrophic assimilation of N
    Hetero_N_assim = (Rh * HGE) / C_Nhetero,
    # Total nitrogen demand
    Ndemand = Auto_N_assim + Hetero_N_assim, # unit should be g N m-2 d-1
    # calculate reach length in m
    reachL = c(((v_mi*10)*w_mi*86.4)/K600_daily_mean),
    Ndemand = if_else(Ndemand < 0, 0.0001, Ndemand), # have demand be essentially zero 
   # Calculate no3 supply 
   NO3_supply = c(((86400*Q_Ls*NO3_mgL_dl)/(w_mi*reachL))/1000),
   # Calculate nh3 supply 
   NH4_supply = c(((86400*Q_Ls*NH4_mgL_dl)/(w_mi*reachL))/1000), ## unit should be g N m-2 d-1
   PO4_supply = c(((86400*Q_Ls*(PO4_ugL_dl/1000))/(w_mi*reachL))/1000) ## unit should be g N m-2 d-1
   )

Plots of Ratios of NH4+- N supply to demand, ratios NO3– N supply to demand

Summary of reach based surface water supply

## [1] 4
## [1] 7
## [1] 38
## # A tibble: 2 × 11
##   catch N_supplym NO3_supplym NO3_supplymin NO3_supplymax NH4_supplym
##   <chr>     <dbl>       <dbl>         <dbl>         <dbl>       <dbl>
## 1 BW         6.88        3.76        0.0526         213.        0.417
## 2 GB        18.5         9.13        0.554           73.6       6.14 
## # ℹ 5 more variables: NH4_supplymin <dbl>, NH4_supplymax <dbl>,
## #   PO4_supplym <dbl>, PO4_supplymin <dbl>, PO4_supplymax <dbl>

Figure 4: Nitrogen supply and demand seasonally at each reach

Plots of Ratios of NH4+- N supply to demand, ratios NO3– N supply to demand colored by DOY

n_s_grid <- ggarrange(
  NO3_supp_plots,
  N_demand_plot1,
  labels = c("b", "c"),
  ncol = 1, nrow = 2,
  label.x = 0.95,  # Move label to the right
  label.y = 1,     # Keep label at the top
  hjust = 1,       # Right-align the label
  vjust = 1        # Top-align the label
)


Ndyn_grid <- ggarrange(
  N_ratio_plots,
  n_s_grid,
  ncol = 2, nrow = 1,
    labels = c("a", ""),
  common.legend = TRUE, 
   widths = c(1.75, 2.25),
  legend = "bottom")

Ndyn_grid

Calculations of flow normalized NO3 and NH4 average daily loads by water year:
Plot for nitrogen cycling by reach location

Figure 5 a-c

Correlations between N uptake and …

    1. biomass
    1. chl-a
    1. GPP
    1. ER
    1. Q
    1. temp
Dataframe for successfull N uptatke

covariat_N <- covariat_datq %>%
  filter(Uadd_ug_L_min < 121) %>%
  filter(success=="yes")%>%
  dplyr::select(site, catch, date, water_year, Chla_ugL_Q, method, Uadd_ug_L_min, 
         Uadd_min_sd, sw_m, sw_sd, Vf_m_min, vf_min_sd, GPP_mean, ER_mean,
         AFDM_ugcm2, Q_m, wtr_m, AFDM_mgg) %>%
    drop_na(Uadd_ug_L_min) 

covariat_NU <- covariat_N %>% distinct()


#hist(covariat_NU$Uadd_ug_L_min,  main = paste(NULL))
#hist(log(covariat_NU$Uadd_ug_L_min+1),  main = paste(NULL))

covariat_NU_B <- covariat_NU%>%
  filter(catch=="BW")

covariat_NU_G <- covariat_NU%>%
  filter(catch=="GB")

U_mod_OM <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(AFDM_mgg) + (1|site), data=covariat_NU_B)
summary(U_mod_OM)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(AFDM_mgg) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 49.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3988 -0.5762  0.1545  0.4827  2.0762 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.6351   0.7969  
##  Residual             0.5798   0.7615  
## Number of obs: 20, groups:  site, 2
## 
## Fixed effects:
##                 Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)       2.6466     0.6030  0.9346   4.389   0.1557  
## scale(AFDM_mgg)  -0.5344     0.1907 17.7468  -2.802   0.0119 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(AFDM_m) -0.087
#hist(residuals(U_mod_OM),main = paste(NULL))

U_mod_bio_r2_vals <- r.squaredGLMM(U_mod_OM)
U_mod_bio_r2_vals
##            R2m      R2c
## [1,] 0.1903139 0.613577

U_mod_OM <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(AFDM_mgg) + (1|site), data=covariat_NU_G)
summary(U_mod_OM)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(AFDM_mgg) + (1 | site)
##    Data: covariat_NU_G
## 
## REML criterion at convergence: 71.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9161 -0.6175 -0.1110  0.3873  2.9432 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.4754   0.6895  
## Number of obs: 33, groups:  site, 2
## 
## Fixed effects:
##                 Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)       2.6164     0.1200 31.0000  21.799   <2e-16 ***
## scale(AFDM_mgg)  -0.3295     0.1219 31.0000  -2.703    0.011 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(AFDM_m) 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_OM),main = paste(NULL))

U_mod_bio_r2_vals <- r.squaredGLMM(U_mod_OM)
U_mod_bio_r2_vals
##            R2m       R2c
## [1,] 0.1859186 0.1859186

d) Uptake and epilithic biomass:

BW:

U_mod_bio <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(AFDM_ugcm2) + (1|site), data=covariat_NU_B)
summary(U_mod_bio)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(AFDM_ugcm2) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 34.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.70519 -0.73386  0.07471  0.87520  1.21436 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.5019   0.7085  
##  Residual             0.6167   0.7853  
## Number of obs: 14, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)         2.5939     0.5743  0.8498   4.517   0.1712  
## scale(AFDM_ugcm2)  -0.4012     0.2211 11.2340  -1.815   0.0964 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(AFDM_2) 0.056
#hist(residuals(U_mod_bio),main = paste(NULL))

U_mod_bio_r2_vals <- r.squaredGLMM(U_mod_bio)
U_mod_bio_conditional_r2 <- round(U_mod_bio_r2_vals[2], 2)

GB:

U_mod_biog <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(AFDM_ugcm2) + (1|site), data=covariat_NU_G)
summary(U_mod_biog)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(AFDM_ugcm2) + (1 | site)
##    Data: covariat_NU_G
## 
## REML criterion at convergence: 56.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55876 -0.69352 -0.08685  0.77110  2.09503 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.7326   0.8559  
## Number of obs: 22, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)        2.46975    0.18249 20.00000  13.534 1.58e-11 ***
## scale(AFDM_ugcm2) -0.01063    0.18678 20.00000  -0.057    0.955    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(AFDM_2) 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_bio),main = paste(NULL))
U_mod_bio_r2_valsg <- r.squaredGLMM(U_mod_bio)
U_mod_bio_conditional_r2g <- round(U_mod_bio_r2_vals[2], 2)

e) Uptake and chl-a:

BW:

U_mod_chla <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(Chla_ugL_Q) + (1|site), data=covariat_NU_B)

summary(U_mod_chla)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(Chla_ugL_Q) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 31.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.20993 -0.01623  0.27721  0.51286  1.33623 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.000    0.0000  
##  Residual             0.664    0.8148  
## Number of obs: 13, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)         2.5657     0.2260 11.0000  11.353 2.05e-07 ***
## scale(Chla_ugL_Q)  -0.5248     0.2352 11.0000  -2.231   0.0474 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(Ch_L_Q) 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_chla),main = paste(NULL))
r.squaredGLMM(U_mod_chla)
##           R2m      R2c
## [1,] 0.293197 0.293197
U_mod_chla_r2_vals <- r.squaredGLMM(U_mod_chla)
U_mod_chla_conditional_r2 <- round(U_mod_chla_r2_vals[2], 2)

GB:

U_mod_chlag <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale((Chla_ugL_Q)) + (1|site), data=covariat_NU_G)

summary(U_mod_chlag)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale((Chla_ugL_Q)) + (1 | site)
##    Data: covariat_NU_G
## 
## REML criterion at convergence: 32.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0086 -0.5137  0.0005  0.7429  1.9490 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.2141   0.4627  
## Number of obs: 22, groups:  site, 2
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)          2.69228    0.09866 20.00000  27.290  < 2e-16 ***
## scale((Chla_ugL_Q)) -0.49926    0.10098 20.00000  -4.944 7.81e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## sc((C_L_Q)) 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_chlag),main = paste(NULL))
r.squaredGLMM(U_mod_chlag)
##            R2m       R2c
## [1,] 0.5379097 0.5379097
U_mod_chla_r2_vals <- r.squaredGLMM(U_mod_chlag)
U_mod_chla_conditional_r2 <- round(U_mod_chla_r2_vals[2], 2)

f) Uptake and GPP:

BW:

U_mod_GPP <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(log(GPP_mean+1)) + (1|site), data=covariat_NU_B)

summary(U_mod_GPP)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(log(GPP_mean + 1)) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 28.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.28637 -0.58032 -0.03618  0.21176  1.89708 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.09332  0.3055  
##  Residual             0.75517  0.8690  
## Number of obs: 11, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               2.13774    0.33992 0.94015   6.289    0.111
## scale(log(GPP_mean + 1))  0.05369    0.27707 8.37351   0.194    0.851
## 
## Correlation of Fixed Effects:
##             (Intr)
## s((GPP_+1)) 0.006
#hist(residuals(U_mod_GPP), main = paste(NULL))

U_mod_GPP_r2_vals <- r.squaredGLMM(U_mod_GPP)
U_mod_GPP_conditional_r2 <- round(U_mod_GPP_r2_vals[2], 2)

GB:

U_mod_GPPg <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(log(GPP_mean+1)) + (1|site), data=covariat_NU_G)

summary(U_mod_GPPg)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(log(GPP_mean + 1)) + (1 | site)
##    Data: covariat_NU_G
## 
## REML criterion at convergence: 9.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0848 -0.2151  0.1633  0.2847  0.8518 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 17.1989  4.1472  
##  Residual              0.1768  0.4205  
## Number of obs: 5, groups:  site, 2
## 
## Fixed effects:
##                          Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)                2.9139     2.9418 0.9007   0.991    0.517  
## scale(log(GPP_mean + 1))   3.2949     0.7635 2.2249   4.315    0.041 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## s((GPP_+1)) 0.045
#hist(residuals(U_mod_GPP), main = paste(NULL))
U_mod_GPP_r2_vals <- r.squaredGLMM(U_mod_GPPg)
U_mod_GPP_conditional_r2 <- round(U_mod_GPP_r2_vals[2], 2)

g) Uptake and |ER|:

BW:

U_mod_ER <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale((ER_mean*-1)) + (1|site), data=covariat_NU_B)

summary(U_mod_ER)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale((ER_mean * -1)) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 25.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7221 -0.4105  0.1035  0.5490  1.3523 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.5886   0.7672  
## Number of obs: 11, groups:  site, 2
## 
## Fixed effects:
##                       Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)             2.1254     0.2313 9.0000   9.188  7.2e-06 ***
## scale((ER_mean * -1))   0.4466     0.2426 9.0000   1.841   0.0988 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## s((ER_*-1)) 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_ER), main = paste(NULL))
U_mod_ER_r2_vals <- r.squaredGLMM(U_mod_ER)
U_mod_ER_conditional_r2 <- round(U_mod_ER_r2_vals[2], 2)

GB:

U_mod_ER <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale((ER_mean*-1)) + (1|site), data=covariat_NU_G)

summary(U_mod_ER)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale((ER_mean * -1)) + (1 | site)
##    Data: covariat_NU_G
## 
## REML criterion at convergence: 7.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8134 -0.6034 -0.1037  0.3973  1.1232 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.09622  0.3102  
##  Residual             0.17679  0.4205  
## Number of obs: 5, groups:  site, 2
## 
## Fixed effects:
##                       Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)             2.3118     0.2904 0.8941   7.960   0.0976 .
## scale((ER_mean * -1))   0.9306     0.2155 2.1657   4.317   0.0430 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## s((ER_*-1)) -0.022
#hist(residuals(U_mod_ER), main = paste(NULL))
U_mod_ER_r2_vals <- r.squaredGLMM(U_mod_ER)
U_mod_ER_conditional_r2 <- round(U_mod_ER_r2_vals[2], 2)

h) Uptake and Q:

BW:

U_mod_Q <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(log(Q_m+1)) + (1|site), data=covariat_NU_B)

summary(U_mod_Q)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(log(Q_m + 1)) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 42.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0014 -0.3472 -0.2333  0.5525  2.2220 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.1599   0.3999  
##  Residual             0.5333   0.7303  
## Number of obs: 18, groups:  site, 2
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)           2.4174     0.3353  0.8893   7.209   0.1075  
## scale(log(Q_m + 1))   0.4530     0.1852 15.8680   2.446   0.0265 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl((Q_+1)) 0.047
#hist(residuals(U_mod_Q), main = paste(NULL))
U_mod_Q_r2_vals <- r.squaredGLMM(U_mod_Q)
U_mod_Q_conditional_r2 <- round(U_mod_Q_r2_vals[2], 2)

GB:

U_mod_Qg <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(log(Q_m)+1) + (1|site), data=covariat_NU_B)

summary(U_mod_Qg)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(log(Q_m) + 1) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 44.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9829 -0.5508 -0.1698  0.4880  2.2779 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.09962  0.3156  
##  Residual             0.61880  0.7866  
## Number of obs: 18, groups:  site, 2
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)           2.3941     0.2946  0.8186   8.127   0.1116  
## scale(log(Q_m) + 1)   0.3448     0.1974 15.9941   1.747   0.0998 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl((Q_)+1) 0.044
#hist(residuals(U_mod_Qg), main = paste(NULL))
U_mod_Q_r2_vals <- r.squaredGLMM(U_mod_Qg)
U_mod_Q_conditional_r2 <- round(U_mod_Q_r2_vals[2], 2)

i) Uptake and temp:

BW:

U_mod_temp <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(wtr_m) + (1|site), data=covariat_NU_B)

summary(U_mod_temp)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(wtr_m) + (1 | site)
##    Data: covariat_NU_B
## 
## REML criterion at convergence: 45
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4477 -0.6101 -0.2310  0.5474  2.1467 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.6829   0.8264  
## Number of obs: 18, groups:  site, 2
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    2.3398     0.1948 16.0000  12.013 2.03e-09 ***
## scale(wtr_m)  -0.2563     0.2004 16.0000  -1.279    0.219    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scal(wtr_m) 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_temp), main = paste(NULL))
r.squaredGLMM(U_mod_temp)
##             R2m        R2c
## [1,] 0.08772759 0.08772759

GB:

U_mod_tempg <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(wtr_m) + (1|site), data=covariat_NU_G)

summary(U_mod_tempg)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(wtr_m) + (1 | site)
##    Data: covariat_NU_G
## 
## REML criterion at convergence: 71.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0146 -0.6432  0.1866  0.6857  2.3690 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.05039  0.2245  
##  Residual             0.63804  0.7988  
## Number of obs: 29, groups:  site, 2
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)    2.5496     0.2175  0.9197  11.720   0.0652 .
## scale(wtr_m)   0.1399     0.1540 26.9996   0.908   0.3718  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scal(wtr_m) 0.010
#hist(residuals(U_mod_temp), main = paste(NULL))
r.squaredGLMM(U_mod_tempg)
##             R2m        R2c
## [1,] 0.02764357 0.09881207

U_mod_NO3 <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(NO3_mgL_dl) + (1|site), data=covariat_datq_BW)

summary(U_mod_NO3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(NO3_mgL_dl) + (1 | site)
##    Data: covariat_datq_BW
## 
## REML criterion at convergence: 183.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.29637 -0.73973  0.05334  0.99156  2.17629 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.0000   0.0000  
##  Residual             0.5996   0.7743  
## Number of obs: 77, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)        2.24906    0.08956 75.00000  25.113   <2e-16 ***
## scale(NO3_mgL_dl) -0.11925    0.07567 75.00000  -1.576    0.119    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(NO3_L_) -0.171
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_NO3), main = paste(NULL))
r.squaredGLMM(U_mod_NO3)
##             R2m        R2c
## [1,] 0.03164349 0.03164349

U_mod_NO3 <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(NO3_mgL_dl) + (1|site), data=covariat_datq_GB%>%filter(NO3_mgL_dl<0.7780))

summary(U_mod_NO3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(NO3_mgL_dl) + (1 | site)
##    Data: covariat_datq_GB %>% filter(NO3_mgL_dl < 0.778)
## 
## REML criterion at convergence: 169.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2861 -0.2404  0.1208  0.4928  2.2231 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.1917   0.4378  
##  Residual             0.8532   0.9237  
## Number of obs: 62, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)        2.52938    0.33122  0.99863   7.637   0.0831 .
## scale(NO3_mgL_dl)  0.05608    0.24008 59.10558   0.234   0.8161  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(NO3_L_) -0.019
#hist(residuals(U_mod_NO3), main = paste(NULL))
r.squaredGLMM(U_mod_NO3)
##               R2m       R2c
## [1,] 0.0007345274 0.1840798

U_mod_NH4 <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(NH4_mgL_dl) + (1|site), data=covariat_datq_BW%>%filter(substrate=="pw"))

summary(U_mod_NH4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(NH4_mgL_dl) + (1 | site)
##    Data: covariat_datq_BW %>% filter(substrate == "pw")
## 
## REML criterion at convergence: 65.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9969 -0.8189 -0.1895  0.7163  2.2768 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.000    0.0000  
##  Residual             0.515    0.7177  
## Number of obs: 29, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)         2.2989     0.1335 27.0000  17.219 4.37e-16 ***
## scale(NH4_mgL_dl)   0.3134     0.1040 27.0000   3.013  0.00556 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(NH4_L_) -0.061
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(U_mod_NH4), main = paste(NULL))
r.squaredGLMM(U_mod_NH4)
##           R2m      R2c
## [1,] 0.244845 0.244845

U_mod_NH4 <- lmer(log(Uadd_ug_L_min+1) ~ 
                      scale(NH4_mgL_dl) + (1|site), data=covariat_datq_GB%>%filter(substrate=="pw"))

summary(U_mod_NH4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Uadd_ug_L_min + 1) ~ scale(NH4_mgL_dl) + (1 | site)
##    Data: covariat_datq_GB %>% filter(substrate == "pw")
## 
## REML criterion at convergence: 94.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9824 -0.4760 -0.0800  0.4451  2.9994 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.06356  0.2521  
##  Residual             0.42799  0.6542  
## Number of obs: 45, groups:  site, 2
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)        2.65248    0.20399  0.91437  13.003   0.0600 .
## scale(NH4_mgL_dl) -0.26077    0.09683 42.35027  -2.693   0.0101 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(NH4_L_) 0.021
#hist(residuals(U_mod_NH4), main = paste(NULL))
r.squaredGLMM(U_mod_NH4)
##            R2m       R2c
## [1,] 0.1394342 0.2507134

Plots of Ut and all linear model covariates by stream

Figure 5: General drivers of either NO3 or NH4 uptake rates (Ut)

Metabolism by water year:

Dataframes:

For metabolism by water year March - October (2021-2024)

Metabolimsm regression analysis :

For nitrogen dynamics by water year March - October (2021-2023)

(3.23-0.79)/3.23
## [1] 0.755418
(13.680257-10.864423)/13.680257
## [1] 0.205832
(0.20-0.08)/0.20
## [1] 0.6
(6.211836-3.511)/6.211836
## [1] 0.4347887

Linear models for metabolism across dry or wet years

lmm_WY_GPP <- lmer(GPP_mean ~ WY_lab + (1|site) + (1|catch), data = metab_WY_lmdf)
summary(lmm_WY_GPP)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GPP_mean ~ WY_lab + (1 | site) + (1 | catch)
##    Data: metab_WY_lmdf
## 
## REML criterion at convergence: 5150.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1777 -0.5184 -0.1201  0.4356  6.6409 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.06227  0.2495  
##  catch    (Intercept) 1.96383  1.4014  
##  Residual             1.67037  1.2924  
## Number of obs: 1531, groups:  site, 4; catch, 2
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     1.77100    1.00001    1.00257   1.771    0.327    
## WY_labnormal   -0.92626    0.08342 1526.89813 -11.104   <2e-16 ***
## WY_labwet      -1.47328    0.08208 1509.10270 -17.950   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.028       
## WY_labwet   -0.032  0.327
#hist(residuals(lmm_WY_GPP), main = paste(NULL))
r.squaredGLMM(lmm_WY_GPP)
##             R2m       R2c
## [1,] 0.09936061 0.5930161
lmm_WY_GPP_BW <- lmer(GPP_mean ~ WY_lab + (1|site), data = metab_WY_lmdf%>%filter(catch=="BW"))
summary(lmm_WY_GPP_BW)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GPP_mean ~ WY_lab + (1 | site)
##    Data: metab_WY_lmdf %>% filter(catch == "BW")
## 
## REML criterion at convergence: 2807.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0656 -0.6096 -0.1146  0.4305  4.6548 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.1986   0.4456  
##  Residual             2.8763   1.6960  
## Number of obs: 719, groups:  site, 2
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    3.3460     0.3287   1.0789   10.18   0.0528 .  
## WY_labnormal  -1.9112     0.1812 715.3341  -10.55   <2e-16 ***
## WY_labwet     -2.6190     0.1460 714.0794  -17.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.129       
## WY_labwet   -0.189  0.285
#hist(residuals(lmm_WY_GPP_BW), main = paste(NULL))
r.squaredGLMM(lmm_WY_GPP_BW)
##            R2m       R2c
## [1,] 0.3257379 0.3692876
lmm_WY_GPP_GB <- lmer(GPP_mean ~ WY_lab + (1|site), data = metab_WY_lmdf%>%filter(catch=="GB"))
summary(lmm_WY_GPP_GB)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GPP_mean ~ WY_lab + (1 | site)
##    Data: metab_WY_lmdf %>% filter(catch == "GB")
## 
## REML criterion at convergence: 172.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.7543 -0.4954 -0.1098 -0.0996 11.0680 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.00000  0.0000  
##  Residual             0.07097  0.2664  
## Number of obs: 812, groups:  site, 2
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    0.20156    0.01352 809.00000  14.903  < 2e-16 ***
## WY_labnormal  -0.17131    0.02153 809.00000  -7.957 5.92e-15 ***
## WY_labwet     -0.11845    0.02445 809.00000  -4.844 1.53e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.628       
## WY_labwet   -0.553  0.347
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(lmm_WY_GPP_GB), main = paste(NULL))
r.squaredGLMM(lmm_WY_GPP_GB)
##           R2m      R2c
## [1,] 0.077609 0.077609
lmm_WY_ER <- lmer((ER_mean*-1) ~ WY_lab + (1|site) + (1|catch), data = metab_WY_lmdf)
summary(lmm_WY_ER)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (ER_mean * -1) ~ WY_lab + (1 | site) + (1 | catch)
##    Data: metab_WY_lmdf
## 
## REML criterion at convergence: 7703.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1922 -0.4857  0.0530  0.5767  4.1823 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 21.673   4.655   
##  catch    (Intercept) 17.285   4.158   
##  Residual              8.827   2.971   
## Number of obs: 1531, groups:  site, 4; catch, 2
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    10.1707     3.7516    1.0008   2.711    0.225    
## WY_labnormal   -0.2992     0.1918 1525.0800  -1.560    0.119    
## WY_labwet      -2.2281     0.1892 1525.2642 -11.779   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.017       
## WY_labwet   -0.019  0.327
#hist(residuals(lmm_WY_ER), main = paste(NULL))
r.squaredGLMM(lmm_WY_ER)
##             R2m       R2c
## [1,] 0.01872544 0.8187402
lmm_WY_ER_BW <- lmer((ER_mean*-1) ~ WY_lab + (1|site), data = metab_WY_lmdf%>%filter(catch=="BW"))
summary(lmm_WY_ER_BW)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (ER_mean * -1) ~ WY_lab + (1 | site)
##    Data: metab_WY_lmdf %>% filter(catch == "BW")
## 
## REML criterion at convergence: 3802.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9074 -0.6550  0.1540  0.6621  3.0292 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 23.64    4.862   
##  Residual             11.49    3.390   
## Number of obs: 719, groups:  site, 2
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   15.0178     3.4430   1.0028   4.362  0.14295    
## WY_labnormal  -1.1740     0.3622 715.0129  -3.241  0.00124 ** 
## WY_labwet     -4.8123     0.2923 715.1599 -16.466  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.025       
## WY_labwet   -0.036  0.284
#hist(residuals(lmm_WY_ER_BW), main = paste(NULL))
r.squaredGLMM(lmm_WY_ER_BW)
##           R2m       R2c
## [1,] 0.119547 0.7119466
lmm_WY_ER_GB <- lmer((ER_mean*-1) ~ WY_lab + (1|site), data =  metab_WY_lmdf%>%filter(catch=="GB"))
summary(lmm_WY_ER_GB)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (ER_mean * -1) ~ WY_lab + (1 | site)
##    Data: metab_WY_lmdf %>% filter(catch == "GB")
## 
## REML criterion at convergence: 3449.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0389 -0.4967 -0.2066  0.5091  5.0202 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 20.045   4.477   
##  Residual              4.039   2.010   
## Number of obs: 812, groups:  site, 2
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    5.3550     3.1678   1.0013   1.690     0.34    
## WY_labnormal   0.6775     0.1644 808.0233   4.121 4.15e-05 ***
## WY_labwet      0.8754     0.1899 808.0543   4.610 4.67e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.022       
## WY_labwet   -0.021  0.370
#hist(residuals(lmm_WY_ER_GB), main = paste(NULL))
r.squaredGLMM(lmm_WY_ER_GB)
##              R2m       R2c
## [1,] 0.006115623 0.8333196

lmm_WY_NO3_Q <- lmer(NO3_supply ~ WY_lab + (1|site) + (1|catch),, data = covariat_WY)
summary(lmm_WY_NO3_Q)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NO3_supply ~ WY_lab + (1 | site) + (1 | catch)
##    Data: covariat_WY
## 
## REML criterion at convergence: 591.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.4708 -0.3057 -0.2494 -0.0662  6.9382 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  24.53    4.953  
##  catch    (Intercept)   0.00    0.000  
##  Residual             826.57   28.750  
## Number of obs: 64, groups:  site, 4; catch, 2
## 
## Fixed effects:
##              Estimate Std. Error     df t value Pr(>|t|)
## (Intercept)     6.288      5.769  4.971   1.090    0.326
## WY_labnormal   -4.460     20.908 59.542  -0.213    0.832
## WY_labwet       6.042      7.780 35.558   0.777    0.442
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.191       
## WY_labwet   -0.610  0.143
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(lmm_WY_NO3_Q), main = paste(NULL))
r.squaredGLMM(lmm_WY_NO3_Q)
##             R2m        R2c
## [1,] 0.01204887 0.04052549
lmm_WY_NH4_Q <- lmer(NH4_supply ~ WY_lab + (1|site) + (1|catch), data = covariat_WY)
summary(lmm_WY_NH4_Q)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NH4_supply ~ WY_lab + (1 | site) + (1 | catch)
##    Data: covariat_WY
## 
## REML criterion at convergence: 352.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8007 -0.4836 -0.2467  0.2983  3.1005 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 31.23    5.588   
##  catch    (Intercept)  0.00    0.000   
##  Residual             44.64    6.681   
## Number of obs: 54, groups:  site, 4; catch, 2
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)    0.4901     3.1800  3.8527   0.154   0.8852  
## WY_labnormal  -3.0047     4.9126 48.1412  -0.612   0.5437  
## WY_labwet      4.9876     1.9978 49.7590   2.497   0.0159 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.097       
## WY_labwet   -0.350  0.153
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(lmm_WY_NH4_Q), main = paste(NULL))
r.squaredGLMM(lmm_WY_NH4_Q)
##             R2m       R2c
## [1,] 0.08667957 0.4625832
lmm_WY_demand_Q <- lmer(Ndemand ~ WY_lab + (1|site) + (1|catch), data = covariat_WY)
summary(lmm_WY_demand_Q)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Ndemand ~ WY_lab + (1 | site) + (1 | catch)
##    Data: covariat_WY
## 
## REML criterion at convergence: -218.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3101 -0.5609 -0.2249  0.4222  3.0156 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  site     (Intercept) 0.0003462 0.01861 
##  catch    (Intercept) 0.0027512 0.05245 
##  Residual             0.0029978 0.05475 
## Number of obs: 81, groups:  site, 4; catch, 2
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   0.08040    0.03940  1.05094   2.040   0.2805    
## WY_labnormal -0.06722    0.03968 75.26618  -1.694   0.0944 .  
## WY_labwet    -0.07771    0.01326 70.13310  -5.858 1.38e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WY_lbn
## WY_labnorml -0.042       
## WY_labwet   -0.169  0.120
#hist(residuals(lmm_WY_demand_Q),main = paste(NULL))
r.squaredGLMM(lmm_WY_demand_Q)
##            R2m       R2c
## [1,] 0.1983456 0.6057187

plots

Figure 6: Average annual differences in GPP, ER, N-demand and N-supply by water year

Values averaged across March to October observations

Extra?

Relationship between Ut and n-demand? not really there.

Correlations btwn GPP and other standing stock estiamtes of productivity

Looking at stoichiometeric ratios of N:P overtime?

N:P

Looks like N:P <16 across all observations

Statisical relationship?

Maybe? GPP and N:P are slightly negatively related. But falls apart when both catchments are modeled independently

mod <- lmer((GPP_mean) ~ scale(N_P_ratio) + (1|catch), data = Stoich_df %>%filter(water_year<2024))
summary(mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (GPP_mean) ~ scale(N_P_ratio) + (1 | catch)
##    Data: Stoich_df %>% filter(water_year < 2024)
## 
## REML criterion at convergence: 163.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4738 -0.4973 -0.1762  0.0920  2.7231 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  catch    (Intercept) 1.121    1.059   
##  Residual             1.551    1.246   
## Number of obs: 49, groups:  catch, 2
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)        0.7355     0.7719  1.0053   0.953   0.5147  
## scale(N_P_ratio)  -0.5650     0.2958 46.0356  -1.910   0.0623 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(N_P_rt) 0.058
#hist(residuals(mod), main = paste(NULL))
r.squaredGLMM(mod)
##             R2m      R2c
## [1,] 0.04253083 0.444168
modb <- lmer((GPP_mean) ~ scale(N_P_ratio) + (1|site), data = Stoich_df %>%filter(water_year<2024 & catch=="BW"))
summary(modb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (GPP_mean) ~ scale(N_P_ratio) + (1 | site)
##    Data: Stoich_df %>% filter(water_year < 2024 & catch == "BW")
## 
## REML criterion at convergence: 110.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3352 -0.5563 -0.3577  0.5690  2.0395 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.2365   0.4863  
##  Residual             2.3756   1.5413  
## Number of obs: 30, groups:  site, 2
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)        1.3925     0.4515  0.9108   3.084    0.219
## scale(N_P_ratio)  -0.7080     0.5118 27.4652  -1.383    0.178
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(N_P_rt) 0.148
#hist(residuals(modb), main = paste(NULL))
r.squaredGLMM(modb)
##             R2m       R2c
## [1,] 0.06103181 0.1460509
modg <- lmer((GPP_mean) ~ scale(N_P_ratio) + (1|site), data = Stoich_df %>%filter(water_year<2024 & catch=="GB"))
summary(modg)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (GPP_mean) ~ scale(N_P_ratio) + (1 | site)
##    Data: Stoich_df %>% filter(water_year < 2024 & catch == "GB")
## 
## REML criterion at convergence: -16.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.6729 -0.5975 -0.4147  0.3845  3.2583 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.00000  0.000   
##  Residual             0.01664  0.129   
## Number of obs: 19, groups:  site, 2
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)       0.07495    0.03047 17.00000   2.459   0.0249 *
## scale(N_P_ratio)  0.02788    0.04747 17.00000   0.587   0.5647  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(N_P_rt) 0.238 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
#hist(residuals(modg), main = paste(NULL))
r.squaredGLMM(modg)
##             R2m        R2c
## [1,] 0.01880352 0.01880352

Statisical relationship?

No

mod1 <- lmer((ER_mean*-1) ~ scale(N_P_ratio) + (1|catch), data = Stoich_df %>%filter(water_year<2024))
summary(mod1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (ER_mean * -1) ~ scale(N_P_ratio) + (1 | catch)
##    Data: Stoich_df %>% filter(water_year < 2024)
## 
## REML criterion at convergence: 253.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8286 -0.6234 -0.2085  0.8830  2.5149 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  catch    (Intercept) 45.06    6.713   
##  Residual             10.06    3.171   
## Number of obs: 49, groups:  catch, 2
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)        8.2966     4.7709  1.0010   1.739    0.332
## scale(N_P_ratio)   0.4253     0.7531 46.0061   0.565    0.575
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(N_P_rt) 0.024
#hist(residuals(mod1), main = paste(NULL))
r.squaredGLMM(mod1)
##             R2m       R2c
## [1,] 0.00121901 0.8177674

Table of mean N:P ratios by site

## # A tibble: 4 × 6
##   site  Tot_N_uMm N_P_ratiom N_P_ratiomin N_P_ratiomax rangeNP
##   <chr>     <dbl>      <dbl>        <dbl>        <dbl>   <dbl>
## 1 BWL       0.559       1.60        0.306         4.11    3.81
## 2 BWU       0.642       3.43        0.223        13.3    13.1 
## 3 GBL       0.805       1.76        0.175         4.91    4.73
## 4 GBU       0.827       1.76        0.332         6.14    5.81

Moddified figure 6?

(3) How does in-stream productivity and nitrogen supply and demand dynamics change within stream catchments and across dry (2021 and 2022) and wet water years (2023)?